function [fars, frrs] = slroc(signals, gtruth, thrs, op)
%SLROC Computes the ROC
%
% [ Syntax ]
%   - [fars, frrs] = slroc(signals, gtruth, thrs, op) 
%
% [ Arguments ]
%   - signals:      the signal values
%   - gtruth:       the groud truths for the signals
%                   It is a logical array of the same size as signals, in 
%                   which 1 indicates that the corresponding signal should 
%                   be accepted, while 0 indicates that it should be rejected.
%   - thrs:         a vector of thresholds at which the rates are computed
%   - op:           the thresholding option
%                   \{:
%                      - high:  accepts if signal >= threshold
%                      - low:   accepts if signal < threshold
%                   \:}
% [ Description ]
%   - [fars, frrs] = slroc(signals, gtruth, thrs, op) computes the
%     false accept rates and false reject rates at specified thresholds
%     for the input signal. 
%
% [ History ]
%   - Created by Dahua Lin on Jun 9th, 2005
%   - Modified by Dahua Lin on May 1st, 2006
%   - Modified by Dahua Lin on Aug 8th, 2006
%       - Base on slhistroc
%   - Modified by Dahua Lin, on Jul 19, 2007
%       - based on new slhistroc
%   

%% parse and verify the input arguments

error(nargchk(4, 4, nargin));

assert(isnumeric(signals), 'sltoolbox:slroc:invalidarg', ...
    'signals must be a numeric array.');

assert(islogical(gtruth), 'sltoolbox:slroc:invalidarg', ...
    'gtruth must be a logical array.');

assert(isequal(size(signals), size(gtruth)), 'sltoolbox:slroc:sizmismatch', ...
    'The sizes of signals and gtruth should be equal.');

assert(isnumeric(thrs) && isvector(thrs) && issorted(thrs), ...
    'sltoolbox:slroc:invalidarg', 'thrs should be a sorted numeric vector.');
cdim = 1 + (size(thrs, 2) > 1);

assert(ischar(op) && (strcmpi(op, 'high') || strcmpi(op, 'low')), ...
    'sltoolbox:slroc:invalidarg', 'op should be either ''high'' or ''low''.');
    

%% compute

if isvector(signals)
    maxv = max(signals);
    minv = min(signals);
else
    maxv = max(signals(:));
    minv = min(signals(:));
end

% case-based dispatch

if thrs(1) <= minv    
    if thrs(end) >= maxv
        edges = thrs;     
        [fars, frrs] = get_frates(signals, gtruth, edges, op);
        
    else
        edges = cat(cdim, thrs, maxv); 
        [fars, frrs] = get_frates(signals, gtruth, edges, op);
        fars = fars(1:end-1);
        frrs = frrs(1:end-1);
        
    end
else
    if thrs(end) >= maxv
        edges = cat(cdim, minv, thrs);
        [fars, frrs] = get_frates(signals, gtruth, edges, op);
        fars = fars(2:end);
        frrs = frrs(2:end);
        
    else
        edges = cat(cdim, minv, thrs, maxv);
        [fars, frrs] = get_frates(signals, gtruth, edges, op);
        fars = fars(2:end-1);
        frrs = frrs(2:end-1);
        
    end
end


%% The core function to compute roc

function [fars, frrs] = get_frates(signals, gtruth, edges, op)

ha = reshape(histc(signals(gtruth), edges), size(edges));
hr = reshape(histc(signals(~gtruth), edges), size(edges));

if isempty(ha)
    ha = zeros(size(edges));
end
if isempty(hr)
    hr = zeros(size(edges));
end

[fars, frrs] = slhistroc(ha, hr, op);

